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Abstract - We establish the connection between a multichannel disordered model -the ID Dirac 
equation with N x N matricial random mass- and a random matrix model corresponding to a de¬ 
formation of the Laguerre ensemble. This allows us to derive exact determinantal representations 
for the density of states and identify its low energy (e —> 0) behaviour p(e) ~ |e|“ _1 . The vanish¬ 
ing of the exponent a for N specific values of the averaged mass over disorder ratio corresponds 
to N phase transitions of topological nature characterised by the change of a quantum number 
(Witten index) which is deduced straightforwardly in the matrix model. 


Since the pioneering work of Dorokhov-Mello-Pereyra- 
Evumar (DMPK) [TJ[2], models of multichannel disordered 
wires have played a prominent role in the theory of dis¬ 
ordered systems as they allow to describe a situation in¬ 
termediate between the strictly one-dimensional (ID) case 
and higher dimensions. With the dimensionality, another 
arucial aspect of disordered systems is the presence or not 
of symmetries, what may strongly affect the localisation 
oroperties. This has led to the classification within or¬ 
thogonal, unitary and symplectic classes, depending on the 
existence of time reversal and spin rotational symmetries, 
according to the Wigner-Dyson classification for random 
matrices, labelled with the Dyson index /3 £ {1, 2, 4}. 
These three classes were later completed by three classes 
with an additional chiral symmetry and four Bogoliubov- 
ie Gennes (BdG) classes for disordered superconductors 
3}[4] (see also ©)■ Several models of multichannel wires 
n these symmetry classes were studied by Brouwer and 
coworkers by extending the DMPK approach [&}|9 . The 
interest for such models has been recently renewed, as 
they can support topologically protected Majorana zero 
modes 10 11 , which are stable against perturbations, 


like a small amount of disorder, what could be used for 
quantum computation [l0 ( 12l. When sufficiently strong, 
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disorder can however drive topological phase transitions 
(quantum phase transitions associated with the change of 
a quantum number of topological nature) 113 15 . 


Although multichannel models can be studied from a 
symmetry viewpoint 16} 17], detailed analysis of specific 
models are also useful, in particular for the determination 
of nonuniversal properties (microscopic parameter depen¬ 
dences) [18]. In this letter, we propose such a detailed 
study of ID models belonging to the chiral classes of disor¬ 
dered systems, and establish the connection with a random 
matrix model defined by the following matrix distribution 


f{Z) =Cn] P (det Z) 


m _1_ / 3(j V _1)/2 


x exp 


--tr { G-^Z + eZ - 1 )} 


( 1 ) 


with matrix argument in the group of N x N Hermitian 
matrices with real (/3 = 1), complex (/? = 2) or quater- 
nionic (/3 = 4) elements and positive eigenvalues. Integra¬ 
tion over this group will be denoted f z>Q DZ f(Z) = 1, 
thus Cjv ,/3 ensures normalisation. We will show that sev¬ 
eral interesting properties of the ID Dirac equation with 
random mass (spectral properties and topological index) 
can be straightforwardly obtained from this matrix dis¬ 
tribution. In the “isotropic case”, when G = g^N, the 
distribution is invariant under unitary transformations. 
Then, Eq. 0 interpolates between the Laguerre distri- 
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bution Ln } o ( Z ) oc (det Z ) 0 ^ 1 exp [ — tr{Z}] with 9 > 0 
(obtained when k —> 0) and the “inverse-Laguerre” distri¬ 
bution In , o { Z ) oc (det 2’)- 6l -/ 3 ( Ar - 1 )- 1 exp [ — tr {Z -1 } ] 
(obtained when k —> oo). The breaking of the invariance 
of /(Z) under unitary transformations, when G is not the 
identity matrix, is a further deformation of the distribu¬ 
tion. Whereas most analytical studies of multichannel dis¬ 
ordered models take as a crucial asumption the isotropy 
among the channels, some of our results will not rely on 
this hypothesis. 

The outline of the letter is as follows : we define our 
model and introduce the scattering problem. Then we ex¬ 
plain how the distribution <[T|) appears. This will provide 
the ground from which localisation and spectral proper¬ 
ties will be recovered directly (in the isotropic case). The 
topological phase transitions, identified from the spectral 
properties, will be characterised through the calculation 
of a topological quantum number, achieved directly from 
a detailed analysis of the matrix distribution (JT|) . 

The model. — The disordered model is the ID Dirac 
equation 'H^(x) = e'S’ix) for 2 N component spinors with 

kL = \a 2 ® 1n9 x + <ti ( 8 ) M(x) , ( 2 ) 

cr, being a Pauli matrix. The Hamiltonian exhibits a chiral 
symmetry a^Hcrs = which puts the problem in one 
of the three chiral classes, depending whether the random 
N x N Hermitian matrix M[x) has real (/3 = 1), complex 
(/3 = 2) or quaternionic (/3 = 4) elements. The mass is 
uncorrelated in space, distributed according to 

P[M(x )] OC e -( 1 / 2 )/ da:tr {(AfG)-MG) t G- 1 (M(x)-/iG)} _ (gj 

The correlations and the mean value ( M(x )) = fiG are 
controlled by the same real symmetric matrix G. The 
dimensionless parameter /i is the averaged mass over dis¬ 
order ratio and will play a central role as it drives the 
topological phase transitions. 

Scattering problem, Riccati matrix and random 
matricial process. — Our starting point is a scattering 
formulation : we consider the Dirac equation on the half 
line with mass M(x) vanishing for x > L in order to set 
a scattering problem. We find convenient to gather the N 
independent solutions of the Dirac equation in a 2TV x N 
“spinor”, 4/ = (</? T ,x T ) T where the two “components” tp 
and x are fw° N x N matrices. In the free region {x > L), 
we can write the spinor as the superposition of incoming 
and outgoing plane waves. For e > 0 we have : 

*(*) = (ilw) e ~ iE{x ~ L) + S ( £ ) « 12 (7i^) e +i£(x “ L) , 

( 4 ) 

where S(e) is the N x N scattering matrix (here charac¬ 
terising the total reflection). Because the chiral symme¬ 
try relates positive and negative energies, we will always 
choose e > 0 (see Ref. 19 for a discussion of the 5-matrix 
symmetry). The study of the strictly ID case (N = 1) 


has emphasized the role of a Riccati variable for provid¬ 
ing the spectral and localisation informations 20-23 . We 


extend here this analysis to the multichannel case and em¬ 
phasize the connection with the scattering problem. We 
define the Riccati matrix as Z e = — ex+A which obeys 
the stochastic differential equation (SDE) : 

d x Z e (x ) = —e 2 — Z E (x) 2 — M(x)Z e (x) — Z e (x)M{x ). (5) 

The initial condition for 4/ [x) is chosen such that the chi¬ 
ral symmetry is preserved : we will choose either i/?(0) = 0 
(which corresponds to M(x) —» — oo on R_) or x(0) = 0 
(M(x) —► +oo on R_).[j] In this case, the Riccati ma¬ 
trix is Hermitian. Matching of Q at the boundary reads 
xiL)tp(L)- 1 = i(l v + 5) (1 jv — S) 1 and allows to express 
the scattering matrix in terms of the Riccati matrix : 


S(e) = [e-iZ e (L)][s + iZ e (L)] 


-l 


( 6 ) 


The SDE © may then be related to a Fokker-Planck 
equation (FPE) describing the matricial random process. 
Without any further assumption (like the isotropy as¬ 
sumption assumed in 3IMI 16]), setting a purely imagi¬ 
nary energy e = i k £ iM, we have obtained (only for /3 = 1 
and 2) that the stationary distribution for the stochas¬ 
tic matricial process is given by Eq. 0 (cf. [39]). When 
N = 1, we check that correspond to the known re¬ 


sults 20 21 


We can make more explicit the relation with the dis¬ 
tribution of the analytically continued scattering ma¬ 
trix S(ik) = {k — Zik)(k + Zjfc) -1 with eigenvalues in 
the interval [—1, +1]. The Jacobian of the transforma¬ 
tion is D Z = D5 (2/c) Ar(1+/3 ( Ar_1 )/ 2 ) det(l - S)~ 2 -^- 1 ), 
thus this shows that the scattering matrix distribution 
is a deformation of the Jacobi distribution : P(S) oc 
det(l + t s)-^-i-/3(^v-i)/2 _ s y-i-p{N-i)/2 exp [- _ 

fctr {G -1 (l + 5 2 )(1 — 5 2 ) -1 } ]. 

Lyapunov spectrum at e = 0. It is well-knwon 
that the introduction of any small amount of disorder, un¬ 
correlated in space, leads to the complete localisation of 
the eigenstates. For weak disorder, |e| g , the localisa¬ 
tion properties of the model 0 are expected to fall into 
the standard universality classes studied in 1 . This does 
not prevent delocalisation to occur at specific points of 
the spectrum as a consequence of some symmetry. The 
Dirac equation may indeed present a delocalisation point, 
exactly at e = 0, as a consequence of the chiral symmetry. 
The energy can therefore be viewed as the chiral-synrmetry 
breaking parameter. 

Localisation can be studied through the concept of Lya¬ 
punov exponents, which measure the exponential growth 
rate of the wave function envelope in different channels. 
For N = 1, the Lyapunov exponent of the model can be 
obtained exactly V e 20. 21||23 . For TV > 1, the determina¬ 
tion of the Lyapunov spectrum is extremely complicated 


1. A more general boundary condition en suring the confinement 


of the particle on M+ is (e 2l0CT2 +0-3) ^(0) = 0 


sin0x(O) = 0. Only 6 = 0 or 7t/2 preserves the chiral symmetry. 
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in general and is only known explicitly at the symmetry 
point (e = 0 ) [g] t hanks to an important simplification^] 
(see also Ref. [16 for a broad perspective on Lyapunov 
spectra in the various symmetry classes). Instead of deal¬ 
ing with the FPE [ 6 }] 8 |, we propose here a more direct 
derivation of the Lyapunov spectrum from the SDE. 

We consider G = gljy- Taking our inspiration from 
the TV = 1 


case 


27 . when e = ik £ 


we introduce 

new variables z n = k , which obey the set of coupled 
SDEs d x ( n = -k sinh 2£ n + gg + & coth(C n - 

Cm) +m ra (x), where fh n (x) are TV independent Gaussian 
white noises of zero mean. The generator of this diffusion 
coincides with the one involved in the FPE of Ref. [§]. For 
jV = 1 . in the absence of the interaction term, we recover 
the SDE of 26 27 . For e = ik = 0, the vanishing of the 


confinment allows for a simple analysis of the dynamics 
as the repulsive forces saturate to ±1. Choosing Ci < 
C 2 < • • • < CJV) one recovers the Lyapunov exponents = 
lim x _ >00 C >n (x)/x of Ref. 6 ] ,[^] directly from the analysis of 
the SDEs : 


In = 


k - - 2 n+ 1 ) 


g for n e {!,••• ,N} . (7) 


Density of states. — The density of states (DoS) of 
multichannel ID models in the chiral and BdG classes was 
studied in several papers when g = 0 [7j[9] : it was shown 
that the low energy DoS exhibits a strong dependence in 
the parity of the channel number TV. Whereas the case of 
even TV leads to a low energy DoS which depends on the 
symmetry index f3, for odd TV the properties of the strictly 
ID Dirac equation with random mass with zero mean were 
recovered, irrespectively of (3 : Dyson singularity of the 
DoS p(e) ~ l/|eIn 3 |e| | . These features were also obtained 
in two BdG classes with broken spin rotational symmetry 
(independently on the parity of TV) [9]. The existence of 
common features for distinct symmetry classes was later 
denoted as “ superuniversality v [29], a concept which was 
used recently for the study of quasi-ID spinless p- wave 
superconducting wires (BdG class D) [l5[|T8] . 

The analysis of the stochastic process Z e (x) allows for a 
direct determination of the spectral properties, what was 
used in the TV = 1 case [20 - 23 . Extending this idea 
when TV > 1, we introduce the characteristic function 
S2 = tr {(Z E (x) + M(x))}, from which the DoS can be 
deduced. The vanishing or the divergence of det Z E (L), 
i.e. of one of its eigenvalues, corresponds to satisfy a 


2. The study of the wave function generally requires to consider 
two coupled random processes controlling the phase and the envelope 
of the wave function. This analysis simplifies either in the high en¬ 
ergy/weak disorder regime when the rapid phase variable decouples 
from the slow envelope variable, or at e = 0 for the Dirac equation, 
because the phase variable is locked and the wave function is only 
controlled by its envelope. 

3. The Lyapunov spectrum 7 ™ = g (N — 2 n + 1) was obtained 

earlier for the matrix T exp [ dx' M(x')] [28| , where T is the 
chronological ordering. Thus, (m) can be understood from the fact 
that the zero energy spinor combines the two independent solutions 
(1, 0) T exp [ fg dx' M(x')\ and (0,1) T exp [ — Jq dx' . 
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Figure 1: DoS p(e) = M'{e) for N = 2 and 3 channels for 
various values of g (mean mass over disorder strength). The 
free DoS is po = N/n. 


second boundary condition for the spinor (detx(-b) = 0 
or det (f{L) = 0) at x = L for one of the eigenmodes. 
As a consequence the integrated DoS per unit length is 
given by = — (1/tt) Im[fi], as for TV = 1 

The mean value {Z e {x)) can be computed from the sta¬ 
tionary distribution f(Z), what is more convenient by 
using that is an analytic function of £ and consider¬ 
ing purely imaginary energy e = ik £ iffi. because f(Z) 
is an equilibrium distribution in this case. Consider¬ 
ing, without loss of generality, G = diag{< 7 i, ■ • • , < 7 tv } 7 we 
deduce the mean value from the normalisation constant 
(\Z E (x) + M(x)\ u ) = gfd\nC Nt p/dgi. 

For the isotropic case, G = g ljv, we write equivalently 

D = (iV^ + fc 511 ^) . (8) 


21 23 


Cn,i 3 can then 
the eigenvalues 


be written as 
of the Riccati 


r dzi "' dZN n i<j k 

„-/i-l-/3(JV-l)/2 


an integral over 
matrix, C/v ,/3 = 

^ where <f>(z) = 

exp[— (z + k 2 3 /z)/(2g)\ 1 which is computed 


using standard technique of random matrix theory 30 
The unitary case (/3 = 2) is the simplest one : we obtain 
the form of a Hankel determinant 


C N , 2 = TV! 2 n k~ N » det [K^+ji-i-^k/g)] , (9) 


where 1 ^ i, j ^ TV. K u (z) is the MacDonald function 


31 


We can deduce exact expressions for the DoS : we 
show in Fig. [T] the DoS for TV = 2 and TV = 3 channels for 
various values of g. The case /3 = 1 is discussed in 39 


Assuming that 0 is also the stationary distribution for 
(3 = 4, we get the Pfaffian 


C N ,4 = N\ 2 n k~ Nfl pf [(j - i) K^ + i + 2 N -i-j (k/ 5 )] , ( 10 ) 
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where 1 ^ i, j ^ 2A. For A = 1, we check that ( |9|10[ ) give 
the known result 20 21 . Similar results were obtained 
for /x = 0 by a slightly different approach in Ref. [9]. 

Using <§ or the equivalent expressions for /3 = 1 and 
/3 = 4, we can analyse in detail the low energy behaviour 
of the DoS, which involves the two leading order terms of 
a k —>■ 0 expansion of Cjv,/ 3 - The DoS presents a power 
law behaviour A /"(e) ~ e a . For /3 = 2, Eq. © shows that 
the next leading order term of the k — > 0 expansion is con¬ 
trolled by an exponent which is a non monotonous function 
of /x, what originates from the expansion of K u (z). The 
exponent is a = 2/x — j3{N — 1) for /x > /3(N — l)/2 and 
vanishes A times below this threshold, presenting a saw 
behaviour (Fig. [2]) (cf. [39]). Each vanishing of the expo¬ 
nent a = 0 corresponds to a critical point where the inte¬ 
grated DoS presents the Dyson singularity J\f(e) ~ 1/ In 2 e 
(the independence in the Dyson index /3 is a sign of su¬ 
peruniversality [29]). When the exponent reaches a local 
maximum (Fig.pl), the behaviour J\f(e) ~ e&\ lne| is found 
(such a behaviour was found in |7|8| for /x = 0 and even A). 




32 


A(/3) = tr{cr 3 exp[-/3-H 2 ]}, 
which can be expressed in terms of the scattering Friedel 
phase <5±(e) = —ilndet5±(£), characterising the scat¬ 
tering problem on the half line x £ R + for (M(x)) = 

±to 0 1wQ : 


A(/3) = 


Z7T 


(ii) 


For a continuous spectrum, A (/?) presents a non trivial /3- 
dependence, however A( 00 ) = [<S+(0) — 5_(0)]/(27t) pro¬ 
vides the information on the number of zero modes. 

As an elementary illustration, we consider first the free 
case (g = 0 ), when the channels are uncoupled and the 


spectrum gapped. Setting M{x) = —00 (i.e. <p(0) = 0) 
and M(x) = tooIat on [0, L], we get the reflection phase 
(for each channel) 8(e )/2 = arg[(fc cosh kL —too sinhftL-|- 
ie sinh/xL)e 17r / 4 j where k = \J trig — e 2 . Letting first 
L —> 00 we deduce that lim^o <$(e) = ndu(mo) + 7 t/ 2 , 
where 0 h(x) is the Heaviside function, leading to A(oo) = 
(A/2) sign(mo)[^] Each channel is characterised by a Wit¬ 
ten index equal to 1 / 2 , which is interpreted as a fermion 
fractionization phenomenon characterising a topologically 
non trivial state [34] . A crucial question is to understand 
how these states are affected by the disorder. Quite re¬ 
markably the zero modes, which are midgap states in the 
absence of disordered, exist in strongly disordered wires 


15 18 with gapless spectrum, as we show below. 


Topological phase transitions. — The spectral den¬ 
sity (bulk information) was related to the mean value 
(Z e ) ; the distribution f(Z) however contains a lot more 
of information. We now determine the topological index 
A(oo) by a direct analysis of ([!]). The relation between 
A(oo) and f(Z) is made more clear by using the supersym¬ 
metry of the Dirac equation : denoting Zf- ( x ) the solution 
of ([5]) for ( M(x )) = ±/xG, we get Z~(x) ll = f — e 2 Z+(x)~ 1 


[21 (equality in law means the same statistical proper¬ 
ties). For e e 1 we deduce (<S + (e) — 8-(e)) /{2 tt) = 
(f/2) En=i (sig n (Zn))- For £ e iK, we get : 


Figure 2: Exponent controlling the low energy integrated DoS 
A f(e) ~ e a as a function of p = mo/g. The dots correspond 
to the critical points where the IDoS is Af(s) ~ l/ln J e. The 
brown squares correspond to the behaviour M{e) ~ e^| lne|. 

Witten index. — We now demonstrate that the van¬ 
ishing of the exponent a is a signature of a sequence of A 
topological phase transitions, accompanied by the change 
of a quantum number of topological origin. From the 
bulk/edge correspondence, we can relate the topological 
number characterising an insulating phase to the number 
of zero modes at the boundary of the system. We con¬ 
sider the Witten index 


1 1 N 

— (8 + {ik) - 8-(ik)) = - ( si g n ( z n - k)) . 


( 12 ) 


n—1 


Below, we will obtain A(oo) from ( 1|12 ). 

For G = glN, the distribution of the Riccati eigenval¬ 
ues for (M(x)) = ±ng 1 N is P±(z lf - ■■ ,z N ) oc n,<j \ z i - 
where (j>{z) = exp[— (z + 

k 2 /z)/2 ] (setting g = 1 for simplicity). In the Coulomb 
gas approach, P± is interpreted as the Gibbs mea¬ 
sure for a ID gas of A “charges” trapped in a con¬ 
fining potential — In <f>(z) and with logarithmic interac¬ 
tions. We show below that the limit k —> 0 involves 
the distribution of eigenvalues of the Laguerre ensemble 
Ajv,e(Ai, • • • , Ajv) = e n,<j I A* — R \ e x ‘/ 2 , 
normalisable for 0 > 0, and the one of the“inverse- 
Laguerre” ensemble 2jv,e(n, ■ • • Tjv) = IL<j l T * “ 

T j\^ Yli T l 6 1 e _1 A 2 ' r i) j i. e . the distribution of 

Tj = 1/Aj’s. We first remark that, for 8 = /x— /3(A— 1)/2 > 
0, the limit k —» 0 of the joint distribution P_(zi, ■ • • , zn) 


corresponds to the Laguerre ensemble. Eq. (12) leads to 


A(oo) = A/2. For p, = p — /3(A — 2n — l)/2 G]0, /?[ with 
n £ {1, • • ■ , A}, the limit k —i 0 is more subtle : we obtain 
that the A charges split into two independent sets, as the 


4. Another connection between topological quantum numbers 
and scattering was discussed in Ref. [33]. 


5. For finite L 1/mo, the midgap state is splitted into two low 
energy states, e± — d=2mo exp[— moL] (for </?(0) = x(L) = 0), which 
is also reflected in the phase shift. Distribution of the ground state 
energy has been obtained in [25[|26| . 
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distribution behaves as 


P-(Z l,"' 


’ ZN h^0 


(Zl Zn\ 

U 2 ’"' ’ k 2 ) 


x Ln — n, /i(-2n+l, 


Z n 

k 2 , 

,Z N ) ■ 


(13) 


(corresponding to the one of Ref. 17]). BdG classes D 


38 


and Dill which support (Majorana) zero modes in ID 
have attracted a lot of interest (see [l7]), hence it would 
be interesting to extend our approach to these classes. 


This distribution describes the “condensation” of n 
charges towards the origin, while the iV — n remaining 
charges are described by a distribution which “freezes” as 
k —> 0. Interestingly, the parameter k drives a phase tran¬ 
sition in the Coulomb gas (the splitting of the charges into 
two groups), which occurs for a finite N, whereas several 
phase transitions were observed up to now in the thermo¬ 
dynamic limit N —» oo in various contexts [35j-[37 . As the 
n smallest charges scale as k 2 —> 0 while the N — n remain¬ 
ing charges scale as fc°, we get from Eq. (12) (Fig. [ 3 ]) : 


A(00) = 


N- 2n 


for ne{0, !,•••, N} . (14) 


Exactly at a critical point g, = j3(N — 2n + l)/2 with n £ 
{!,■•• ,1V}, the charge n must be considered separately : 


,z N ) 


1 


fc -+0 jfc 2 (™~ 1 ) 


T I 21 

-i-n-l.fl I — 


>*( 


k 2 ' 


Zn -1 

k 2 


JL e - ( fe /2 )[*„/fc+*/*»] c N _ n ^ Zn+1 , . .., ZN ). 


(15) 


The position of the charge n now scales as k —> 0, however 
its distribution is symmetric with respect to z n = k, which 
shows that (sign(z n — k)) = 0, thus it does not contribute 
to (12). As a result A(oo) = (N — 2n + l)/2 (Fig. [ 3 ]). A 
fine tuning of the disorder (or the mass) can thus change 
the parity of the number of zero modes. 



Figure 4: Phase diagram of the model with N channels : the 
half plane (disorder against averaged mass) is splitted in N + 
1 sectors characterised by the value of the topological number 
2A(oo). The N lines correspond to g = mo/g = (/3/2)(N — 
2 n + 1), with n £ {1, • - • ,1V}. On each transition line, the 
index takes the intermediate value. 


Although most of the discussion has concerned the 
isotropic case, the main conclusions can be extended to 
the non-isotropic case, G = diag((q, • • • , g^r). The deter- 
minantal representations ( 9|10 ) can be generalised. Using 
these generalisations and the distribution ([!]), we were able 
to show that both the exponent a and the Witten index 
present the same dependence with /i as in the isotropic 
case (Figs. [ 2 ] and [3]) . The Lyapunov analysis is much more 
difficult to extend. Using a formalism based on the QR de¬ 
composition, we have obtained analytical expressions for 
the Lyapunov exponents (at e = 0) in the non-isotropic 
case for N = 2 channels : 


7i = 2 


9i92 
gi + 32 


P 2 (9i92T 


gi +92 


2 p 

9 i 


2 fi 

92 


(16) 


2 

1 


< 

-1 
-2 

"-3 -2-10 1 2 3 

/j ( mass/disorder) 

Figure 3: Witten index A(oo) (number of zero modes is 
2|A(oo)|^. Dashed lines are the Lyapunov exponents. 

This analysis can be confronted to the Lyapunov 
spetrum 0 plotted in Fig. [3]: this suggests that the cre¬ 
ation or the annihilation of a pair of zero modes at a topo¬ 
logical transition is made possible by the delocalisation in 
one of the N channels. 

Conclusion. — In this article we have shown a con¬ 
nection between multichannel ID disordered models with 
chiral symmetry and the random matrix model defined by 
Eq. ([I]), what has allowed us to derive several exact re¬ 
sults straightforwardly (DoS and topological index count¬ 
ing the number of chiral zero modes). This analysis can 
be summarized on the phase diagram plotted in Fig. [4] 
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Where B (z\ a, b ) is the incomplete Beta-function, and Fi is 
the Appell Hypergeometric function with two arguments: 


Fi (a, b , b', c; x, y) = ^ ^ - 


ra^O n^O 


(&)m+n(^)m(^ )r 

(c) m+ „m!n! 


-x"V 


(17) 

with (a) n = T(a + n)/T(a). The second Lyapunov expo¬ 
nent follows from the sum rule 7n = 9 tr{G} (valid 
at e = 0), i.e. 72 = p (3i + 32 ) — 7 i- For N > 2, we have 
performed a numerical analysis. Our results also exhibit 
the vanishing of one Lyapunov at each phase transition, 
as expected (Fig. [5|. Details will be published elsewhere. 

Although we are providing here the first exact results 
for a non-isotropic multichannel disorder, to the best 
of our knowledge, we stress that we have not consid¬ 
ered here the most general situation, which would in¬ 
volve the distribution of disorder of the form P[M(x)\ oc 
exp [ — f da; tr (x)B~ 1 M(x)} ]. This remains a 

challenging problem. 
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Figure 5: Witten index A(oo) and Lyapunov spectra with 
anisotropic disorder : for two channels (gi/g 2 = 3/ and three 
channels (gi/g 2 = 52/53 = 3/. The squares are the values 
obtained numerically. 
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Supplementary material 

Aurelien Grabsch and Christophe Texier 

Abstract — We present some technical details : ( i ) proof 
of the stationary matrix distribution (main result of the let¬ 
ter, Eq. 1). ( ii ) Derivation of the exponent of the low energy 
DoS (the orthogonal and symplectic cases are also discussed). 
Finally we discuss the relation between the Lyapunov analysis 
and the conductance. 

A. The Fokker-Planck equation for the matricial 
process. — Let us write the Fokker-Planck equation as¬ 
sociated to the SDE satisfied by the Riccati matrix : 

Z' = -Z 2 - E — fj,GZ - llZG - ZM - MZ, (18) 

where E = e 2 and M[x) is a Gaussian white noise dis¬ 
tributed according to : 


P[M(x)] oc exp 


dx tr j.M(x)^G 1 M(x)| 


<» "("- 1} + iv, 


g a if a = 6 , 
a ab - 1 iia^b. 


k<l 


where 

Bmn,kl(Z) =- 


One can now write the adjoint generator of the diffusion 
as : 


# = ^ AA— [ z2 + e + VGZ + [iZG] 
OZj n - 

m^cn 

i 

+ 2 


E/ E/E/ f)Z E rnn,kl(Z) gg ^B abk l{Z). 


< K, < h dZ " 
m^n k^.1 a^b 

To take into account the symmetry of Z , introduce the 


differential operator M82 


d_ 

dZ 


1 + Sn 


d 


t 2 dZ mn 
and also a modified matrix of variances : 

: _ 2 -S ab 

Gab — 2 ® cib ’ 


(25) 


(26) 


The generator can then be expressed in the compact form : 


, (19) 


with M a Hermitian NxN matrix with real (/? = 1), com¬ 
plex (j3 = 2) or quaternionic (/3 = 4) elements. The Her¬ 
mitian matrix G controls the correlations. The three cases 
correspond respectively to the classes BDI (chiral orthog¬ 
onal), AIII (chiral unitary) and CII (chiral symplectic) of 
the classification of disordered systems (cf. [EM08] ). One 
can assume without loss of generality that G is diagonal 
G = diag(gi, • • ■ , g^). This implies that the entries of M 
are independent Gaussian white noises. The number of 
independent components of M (and also of Z) being 



d v 


d rr 

Id \ T 1 

11 

2 tr < 

~dZ Z 

a o 

dz z + 

1 

E 

IN 

|co 



+ tr ^ iyZ 2, + E + g,GZ + g,ZG) ^ , 

where [A o B] mn = A mn B mn is the Hadamard product. 

A.2. Case f3 = 2 (class AIII). In this case, one 
can treat the non diagonal entries Z mn and Z nm ^inde- 
pendently. But now, the non diagonal entries of M are 
complex. Denote : 


yjomn(.(mn T l(mn)i 


(27) 


( 20 ) 


where a is given by (22), ( mn and £ mn are Gaussian white 


one has to distinguish the cases /3 = 1 and /3 = 2 in the 
derivation of the Fokker-Planck equation. 

A.l. Case (3 = 1 (class BDI). In this case the ma¬ 
trices M and Z are real and symmetric, so one has to 
consider only their upper-triangular part. Denote 

M ab = yCGbCab, ( 21 ) 

where ( ab is a Gaussian white noise of variance 1 and 


( 22 ) 


9a+9b 

Writing the SDE in components 

Z' mn = [- -Z 2 — E — g,GZ — nZG] mn + J2 

(23) 


(24) 


X (Z mk 8 n l -(- ZmiSnk T Zi n 5 k m T ■ 


noises of variance 1 , with ( mn = £ nm and £ mn = -£„ m . 
The SDE gives, in components : 

Z' mn =[-Z 2 -E-nGZ-nZG\ mn 

+ 'y ] B mn ^i{Z) £ki + E Cmn,kl (Z)tku (28) 

k^l k<l 

with B(Z) given by p4| ), and 

Cmn,kl(Z ) = — y^Oki (29) 

X (Zmk^nl ZmiSnk T Zi n 6 krn Z kn 8im) • 

The generator is now : 

= E J- A + E + f*CZ + l*ZG] mn 

m,n mn 

+oEEEhf Bmn,kl (Z) y B a b : kl (^) 

z ^, , 'Jzjmn Uzjab 

m,n k ^.1 a,b 

+ 2 E E E g^c mn M z )g^-c ab Mz)- 

m,n k<l a,b 
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Denoting simply 


d_ 

dZ 


d 

d. Z™ 


(30) 


We recognize the MacDonald function, hence 
Cn, 2 = N\ 2 N k~ Ntl det [K^+N^^k/ 


(37) 


and again a defined by ( 261 , the generator takes the form : 


# = tr 



-1 T 


G O 


§z ]z 


d 7 t 
dZ Z 


• Low energy analsyis .— The analysis of the k —> 0 be¬ 
haviour of the normalisation constant relies on the prop¬ 
erties of the MacDonald functions : 

, » / k~^(l + A„k 2 + o(k 2 )) for \v\ ^ 1, 

I ' W ~\ k~^{l + B u k 2 ^+o(k 2 ^)) for \v\ < 1 . 


+ 


<9z 


[Z 2 + £7 + ^GZ + /iZG] 


Introduce 


I 


v = n — TV + n e]0,1[, 


A.3. Stationary solution. The well-known form for 


the stationary solution in the case TV = 1 BCGL90 
CTT11] has led us to propose the distribution 


/(Z)=C^det(Z) 


1-j9(JV-1)/2 


x exp 


--trjG-^Z-^Z- 1 )} 


(31) 


Using the expression of the generators given above and 
standard formulae for derivatives of traces and determi¬ 
nants PP12 , we have verified after lengthy calculations 


that (31) is a solution of the stationary Fokker-Planck 
equation : 

#7 = 0. (32) 

B. DoS in the isotropic case. — When the corre¬ 
lation matrix G is proportional to unity, G = gl ^, the 
normalisation constant reduces to a multiple integral over 
the TV eigenvalues of the Riccati matrix : 


(38) 

(39) 

with n £ {1,..., TV}. In the expansion of the determinant, 
many terms give a contribution to the two leading orders. 
To get the power of k, one only needs to exhibit one of 
these terms. For example, the product of all diagonal 
terms gives a contribution to the leading order , giving a 
behaviour k * 1 (the exact expression of g is not important). 
To get the next leading term, one has to distinguish two 
cases depending on the parity of n. 

• Odd n .— On the diagonal, the MacDonald function 
with the smallest index is 17 , which gives a contribu¬ 
tion fc _ "( 1 + B^k 2 "). This gives the expansion Cn ,0 ~ 
k v ( 1 + Ak 2 "), from which one gets 


= 2 v = 27 — TV + n) . 


(40) 


• Even n .— Now, the MacDonald function with the 
smallest index is K v _\. Proceeding as before, it gives 
C N ,p ~ k r >(l + Bk 2 < 1 - v '>), so : 


Cn,p — 


dzi ■ ■ ■ dz N JJ — Zj\P U <j>{z {), (33) 


a = 2(1 — v) = 2(TV - n + 1 - At) • 


(41) 


*< j 


where = z ** 1 1 ^ 2 exp[— (z + k 2 /z)/( 2 g)]. 

Such integrals may be written under the form of deter¬ 


minants thanks to standard technique M04 


The two behaviours ( 40|41 ) correspond to the saw be¬ 
haviour shown in Fig. 2 of the letter. 

B.2. Case /3 = 1 (class BDI). We wish to calculate 
the matrix integral (33) for (3 = 1. The case of even and 


The low energy behaviour of the DoS can be obtained 
from the two leading orders of a k —» 0 expansion of C/v ,/3 : 


odd TV must be treated separately. 

• Even TV.— We use a standard result due to de 


Bruijn B56 which allows to express the normalisation in 


Civ, / 3 ~fc ,) (l + Afc a +o(A; a )). 


(34) 


Using Eq. ( 6 ) of the letter, one gets the expansion of the 
characteristic function : 


S2 = —N p,g — r/g — aAgk a + o(k a ). 


(35) 


Then, analytic continuation to k = —ie yields Af(e) ~ e a . 
We now discuss how the exponent a can be determined. 

B.l. Case (3 = 2 (class AIII). The unitary case is 


terms of a Pfaflian 

Cn,i = TV! pf(A) (42) 

where A is the TV x TV antisymmetric matrix with elements 

POO 

Aij = / dz!dz 2 signal—2 2 ) 7 _1 ~2 _ V(~i)<K~2)- (43) 


Odd TV.— Eq. (42) holds provided one adds a line and 


the most simple one : the integral ( 33 ) may be rewritten a column to the matrix A , given by : 


under the form of a determinant Cn ,2 = IV! det#), where 
the TV x TV matrix A has elements 


A it at +1 = —.Ajv+i,* = / dzz 1 1 <f>(z ) (44) 

Jo 


Ai,j — 


dz <(>(z 


J+j-2 


(36) 


= 2k^ i -—K^_ ti _ i (k). 


(45) 
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• Low energy analysis. — The leading order as k —> 0 can 
be obtained directly from (431 by rescaling either z\ or 2 2 
by a factor k 2 in order to obtain convergent integrals in 
the limit k —)• 0. To get the next order, we rescale the 
variables as z\ = kx and z 2 = ky. This gives : 


B.3. Case f3 = 4 (class CII). Although we have not 
explicitly proven that pT] is true for the case of matrix 
with quartionic elements (/3 = 4), we can compute the 
multiple integral ( [33]) in this case. It can be related to a 
Pfaffian Cn,a = TV! pf(A) of the 2 TV x 2 TV matrix 


Aid = 


(46) 


where we have introduced : 


/>00 

Aij = (,j -i) dz cj>(z) z z+j ~ 3 . (57) 

Jo 


B a ,p(k) = dx dy sign(z - y) (47) 

Jo Jo 

x x a ~ 1 y^~ 1 e~^^ x+y+1 ^ x+1 / y ) 
From this expression, it can be shown that 

(a + (3)B at p(k) + kd k B a ^(k) = - kB a ^i t p(k) 

- kBaf-^k), (48) 

(a - /3)B at p(k) + kd k B at p(k) = - kB a ^ +l (k) 

-kB a _ hP (k) ~ 4K a+ p(2k) . (49) 

These relations can be used to obtain the behaviour of 
B a ,p(k) as k —► 0 : 

B a< p(k) ~ 1 + 0(fc"(“’«)) , (50) 

where 

£(a,P) = max(|a + /3|,|a- (3\) , (51) 


Some algebra gives 

Cjv , 4 = N\ 2 n k~ Nl * (58) 

X pf [O' - i) Ku.+l+ 2 N-i-j(k/g)\ KijK2A r ■ 

The formula obviously give the known result for TV = 1, 
as it should. For TV = 2 we check that ( [58] ) leads to the 
low energy behaviour J\f(e) ~ e 4 ln(l/e) and the expected 
power law (Fig. [ 6 ]). 



and 

V(a,f3) = 2min(l, \a\, |/3|) . (52) 

For odd TV we must add one line and one column to the 
matrix A defined for even TV, as explained above. Both 
cases can then be treated on the same footing. 

Denote 

u = fj,~ (TV-2n- l)/2 e]0,l[ (53) 

with n £ {1,..., TV}. It is easier to compute det(A) by 
decomposing the matrix in blocks (for both odd and even 
TV): 



A = 

( A n 
{-Bin 

Cn,u J 

where we 

introduced 



- A n of size n x n, 



BN,n 

of size n x 

( N-n ), 


~ C N , n 

of size (TV - 

1 

X 

-n). 


One can then obtain the expansion of det(A) as k —> 0 
using the formula for the determinant of a block matrix : 

det(A) ~ det(CW,n) det(A n + B N ^ n C~^ n B N,n) ■ ( 55 ) 

And finally the behaviour of Cjv.i follows from pf(A) = 
\J det(A) ~ fc <T (l + 0 (fc“)), where the exponent controlling 
the next leading order term is 

a = 2 min(i/, 1 — v) . (56) 

This corresponds to the saw behaviour described in the 
body of the letter (curve of Fig. 2 of the letter up to a 
rescaling, see below). 


ei g 

Figure 6: DoS for the symplectic case for N = 2 channels. 

B. f. Low energy DoS - Summary. The previous anal¬ 
ysis has demonstrated that the exponent a controlling the 
low energy intergated DoS, N(e) ~ e a , presents the saw 
behaviour as a function of fi. The dependence in the Dyson 
index (3 can be simply written under the form 

a p(D) = f “2 Qm) • (59) 

In the case /3 = 4 (chiral symplectic case), we have checked 
the relation with ( [58] ) for TV = 2. 

C. Lyapunov spectrum and conductance at a 
topological phase transition. — Away from the TV 
topological phase transitions, all the Lyapunov exponents 
are finite. In this case we expect conventional localisa¬ 
tion properties, reflected for example in the conductance 
Q characterising the transmission through a slice of ran¬ 
dom medium (i.e. the random mass is non zero on [0, L\ 
and vanishes on ] — oo, 0] U [L, +oo[). The conductance is 
expected to present the generic behaviour for ID or quasi- 
ID systems : exponential decay with L. The logarithm of 
the conductance is expected to be self averaging, thus the 
distribution of In Q is expected to be Gaussian, centered 
on (lnt?) ~ — 2 L/£i oc , where £i oc is the localisation length. 

Precisely at a topological phase transition, when the ex¬ 
ponent a vanishes, the model presents atypical localisation 
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properties. For e = 0, the conductance can be expressed in 
terms of the variables as Q = J2 n coslD 2 (, n (L) introduced 
in the paper, where each variable C n (x) behaves asymp¬ 
totically as a free Brownian motion with a drift given by 
the corresponding Lyapunov exponent. When one Lya¬ 
punov exponent vanishes, say 7 „ = 0 , the conductance is 
dominated by the contribution of the corresponding vari¬ 
able thus argcosli[l/ \[Q\ is a Gaussian variable with 
variance gL , as in ID [SCFG 99 (this is a new manifesta¬ 
tion of superuniversality at the topological phase transi¬ 
tion). For L —> 00, the moments of the conductance then 
behave as (Q n ) ~ 2 2n /(n^/2ngL), indicating very large 
fluctuations, with non self averaging In Q when L —> 00 : 
(lnC?) ~ — yj8gL/n and var(ln Q ) ~ 4(1 — 2/n)gL (this was 
obtained in [BFGM00 BMFOl] for /i = 0 and N odd ; see 
also [RM14| for the energy dependence in the case N = 1). 
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